'''
Created on 23.2.2011

@author: Martin Vegi Kysel
@param x,y: two indepentent lists of values
@summary: Implements the Mann Whitney non-paired U test
@return: [p,z,rho,U]
'''

import math
from PyStats.WilcoxonFunction import CumulativeNormalDistribution

criticalUTable05 = [
                        [-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1],
                        [-1,-1,-1,-1,-1,-1,-1,0,0,0,0,1,1,1,1,1,2,2,2,2],
                        [-1,-1,-1,-1,0,1,1,2,2,3,3,4,4,5,5,6,6,7,7,8],
                        [-1,-1,-1,0,1,2,3,4,4,5,6,7,8,9,10,11,11,12,13,13],
                        [-1,-1,0,1,2,3,5,6,7,8,9,11,12,13,14,15,17,18,19,20],
                        [-1,-1,1,2,3,5,6,8,10,11,13,14,16,17,19,21,22,24,25,27],
                        [-1,-1,1,3,5,6,8,10,12,14,16,18,20,22,24,26,28,30,32,34],
                        [-1,0,2,4,6,8,10,13,15,17,19,22,24,26,29,31,34,36,38,41],
                        [-1,0,2,4,7,10,12,15,17,21,23,26,28,31,34,37,39,42,45,48],
                        [-1,0,3,5,8,11,14,17,20,23,26,29,33,36,39,42,45,48,52,55],
                        [-1,0,3,6,9,13,16,19,23,26,30,33,37,40,44,47,51,55,58,62],
                        [-1,1,4,7,11,14,18,22,26,29,33,37,41,45,49,53,57,61,65,69],
                        [-1,1,4,8,12,16,20,24,28,33,37,41,45,50,54,59,63,67,72,76],
                        [-1,1,5,9,13,17,22,26,31,36,40,45,50,55,59,64,67,74,78,83],
                        [-1,1,5,10,14,19,24,29,34,39,44,49,54,59,64,70,75,80,85,90],
                        [-1,1,6,11,15,21,26,31,37,42,47,53,59,64,70,75,81,86,92,98],
                        [-1,2,6,11,17,22,28,34,39,45,51,57,63,67,75,81,87,93,99,105],
                        [-1,2,7,12,18,24,30,36,42,48,55,61,67,74,80,86,93,99,106,112],
                        [-1,2,7,13,19,25,32,38,45,52,58,65,72,78,85,92,99,106,113,119],
                        [-1,2,8,14,20,27,34,41,48,55,62,69,76,83,90,98,105,112,119,127]
                    ]
    
criticalUTable01 =[
                        [-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1],
                        [-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,0,0],
                        [-1,-1,-1,-1,-1,-1,-1,-1,0,0,0,1,1,1,2,2,2,2,3,3],
                        [-1,-1,-1,0,0,1,1,2,2,3,3,4,5,5,6,6,7,8],
                        [-1,-1,-1,0,1,1,2,3,4,5,6,7,7,8,9,10,11,12,13],
                        [-1,-1,-1,0,1,2,3,4,5,6,7,9,10,11,12,13,15,16,17,18],
                        [-1,-1,-1,0,1,3,4,6,7,9,10,12,13,15,16,18,19,21,22,24],
                        [-1,-1,-1,1,2,4,6,7,9,11,13,15,17,18,20,22,24,26,28,30],
                        [-1,-1,0,1,3,5,7,9,11,13,16,18,20,22,24,27,29,31,33,36],
                        [-1,-1,0,2,4,6,9,11,13,16,18,21,24,26,29,31,34,37,39,42],
                        [-1,-1,0,2,5,7,10,13,16,18,21,24,27,30,33,36,39,42,45,46],
                        [-1,-1,1,3,6,9,12,15,18,21,24,27,31,34,37,41,44,47,51,54],
                        [-1,-1,1,3,7,10,13,17,20,24,27,31,34,38,42,45,49,53,56,60],
                        [-1,-1,1,4,7,11,15,18,22,26,30,34,38,42,46,50,54,58,63,67],
                        [-1,-1,2,5,8,12,16,20,24,29,33,37,42,46,51,55,60,64,69,73],
                        [-1,-1,2,5,9,13,18,22,27,31,36,41,45,50,55,60,65,70,74,79],
                        [-1,-1,2,6,10,15,19,24,29,34,39,44,49,54,60,65,70,75,81,86],
                        [-1,-1,2,6,11,16,21,26,31,37,42,47,53,58,64,70,75,81,87,92],
                        [-1,0,3,7,12,17,22,28,33,39,45,51,56,63,69,74,81,87,93,99],
                        [-1,0,3,8,13,18,24,30,36,42,46,54,60,67,73,79,86,92,99,105]                   
                    ] 
    



def tiedRanks(val):
# computes the ranks for a given sample array
# positive and negative samples are combined into one
# the original position in the array remains and is replaced with the rank

    # sort the input array
    sortedVals = sorted(val)
    
    
    # create the output array and fill it with original values
    normalizedVals = [x for x in range(len(val))]


    
    current = sortedVals[0] # start the iteration with the first element
    clusterSize = 0 #the amount of equivalent values ~(7,7,7)
    
    
    for i in range(len(sortedVals)):
        # iterate over all sorted values
        
        # if the current value in the array is the same as the previous
        # increase the clustersize and continue to next value
        if current==sortedVals[i]:
            clusterSize=clusterSize+1
        else: 
            # if the current value is different from the previous
            # compute the mean and exchange the values in the normalized array
            if clusterSize>1:
                sum = 0.0
                
                for j in range(clusterSize):
                    sum=sum+i-j
                    
                # compute the mean I value
                normalizedI = sum/clusterSize
                
                # exchange all values in the array for mean I
                for j in range(len(val)):
                    if(val[j]==current):
                        normalizedVals[j]=normalizedI 
                
                # the last value was different
                clusterSize=1
                
                #choose the next value
                current = sortedVals[i]
            
            else:
                # the current value was found only once in the array
                # compute no mean and exchange for I
                for j in range(len(val)):
                    if(val[j]==current):
                        normalizedVals[j]=i 
                clusterSize=1
                current = sortedVals[i]
                
        # because the logic of this function writes the values AFTER the next is read
        # if the current value is the last, repeat the last steps        
        if i==len(sortedVals)-1:
            if clusterSize>1:
                sum = 0.0
                for j in range(clusterSize):
                    sum=sum+i-j+1
                normalizedI = sum/clusterSize
                for j in range(len(val)):
                    if(val[j]==current):
                        normalizedVals[j]=normalizedI 
            else:
                for j in range(len(val)):
                    if(val[j]==current):
                        normalizedVals[j]=i+1
        
    # return the ranked array    
    return normalizedVals
    
    
    
def MannWhitney(x,y):
# implements the Mann-whitney U test for unpaired samples
# x,y - two independent samples, can be of varying sizes
    
    returnArray = [0,0,0,0]
    
    
    # n1,n2 - sample size, n1 is always the bigger one
    if x>y:
        n1=len(x)
        n2=len(y)
    else:
        n1=len(y)
        n2=len(x)
    print "length of n1= ",n1
    print "length of n2= ",n2

    
    if n1<=20 and n2<=20 and criticalUTable05[n1-1][n2-1]==-1:
        print "sample size too small!"
        return [-1,-1,-1,-1]
    
    # append all samples into one array, start with the bigger array
    allValues = []
    
    if x>y:
        for current in x:
            allValues.append(current)
        for current in y:
            allValues.append(current)
    else:
        for current in y:
            allValues.append(current)
        for current in x:
            allValues.append(current)
    
    # compute the ranks for the samples
    ranks = tiedRanks(allValues)
    
    # rank summaries
    n1Sum = 0;
    n2Sum = 0;
    
    # summarize the ranks for n1 and n2
    for i in range(n1):
        n1Sum=n1Sum+ranks[i]
    for i in range(n2): # start from the last sample of n1
        n2Sum=n2Sum+ranks[i+n1]
    
    
    # compute the U values
    U1 = n1Sum - (n1*(n1+1))/2
    U2 = n2Sum - (n2*(n2+1))/2
    print "U1 ",U1, "U2", U2
    
    # the smaller U is the outcome
    finalU = 0;
    if U1>U2:
        finalU = U2
    else:
        finalU= U1
    
    
    print "Mann-Whitney U-Test. U=",finalU
    
    
    # compute rho
    rho = 1- finalU/(n1*n2)
    print "Rho Value for U(",finalU,") is ",rho
    
    
            
    
    p = 0.0
    
    if n1<=20 and n2<=20:
        # sample size too small for z-approximation
        if criticalUTable05[n1-1][n2-1]==-1:
            print "Sample size too small at 0.05"
        else:
            if finalU<=criticalUTable05[n1-1][n2-1]:
                print "Hypothesis H0 can be rejected at .05",finalU,"<=",criticalUTable05[n1-1][n2-1]
                p = 0.05
            else:
                print "Hypothesis H0 cannot be rejected at .05",finalU,">=",criticalUTable05[n1-1][n2-1]
                #p = 0.0
                
        if criticalUTable01[n1-1][n2-1]==-1:
            print "Sample size too small at 0.01"
        else:
            if finalU<=criticalUTable01[n1-1][n2-1]:
                p = 0.01
                print "Hypothesis H0 can be rejected at .01",finalU,"<=",criticalUTable01[n1-1][n2-1]
            else:
                print "Hypothesis H0 cannot be rejected at .01",finalU,">=",criticalUTable01[n1-1][n2-1]
        
        print "WARNING: following P value must not be exact. Refer to the above mentioned p threshold."
    
        
    upperPartZ = finalU - (n1*n2)/2
    lowerPartZ = math.sqrt((n1*n2*(n1+n2+1))/12)

    z = upperPartZ/lowerPartZ
    print "z score =", z
    
    
    # use the function from Wilcoxon
    print "p (uncorrected, one-tailed)=", CumulativeNormalDistribution(z)
    print "p (uncorrected)=", CumulativeNormalDistribution(z)*2.0
    
    
    if n1>20 or n2>20:
        p = CumulativeNormalDistribution(z)*2.0
    
    
    #print returnValue
    returnArray[0] = p
    returnArray[1] = z
    returnArray[2] = rho
    returnArray[3] = finalU
    
    return returnArray

# to start this, use something as follows:    
#x = [1, 2, 3, 4, 5, 6, 7, 8,9,10,11,12,13,14,15,16,17,18,19,20,21]
#y = [15, 18, 22, 16, 18, 19, 17, 10,9,10,11,12,13,14,15,16,17,18,19,20,21]
#MannWhitney(x, y)
